function low_rank_approximation()
% LOW_RANK_APPROXIMATION 低秩近似算法
% 
% 实现和比较各种低秩近似方法

fprintf('=== 低秩近似算法 ===\n\n');

%% 1. SVD最优低秩近似
fprintf('1. SVD最优低秩近似\n');
fprintf('目标：理解Eckart-Young定理\n\n');

% 构造测试矩阵
m = 20; n = 15; true_rank = 5;
rng(42);
U_true = orth(randn(m, true_rank));
V_true = orth(randn(n, true_rank));
S_true = diag([10, 8, 6, 4, 2]);
A = U_true * S_true * V_true' + 0.1 * randn(m, n);

fprintf('测试矩阵: %d×%d，真实秩≈%d\n', m, n, true_rank);

% SVD分解
[U, S, V] = svd(A);
singular_values = diag(S);
fprintf('前8个奇异值: [%.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f]\n', ...
        singular_values(1:8));

% 不同秩的最优近似
ranks = 1:8;
fprintf('\nSVD最优近似分析:\n');
fprintf('秩    Frobenius误差    2-范数误差    相对误差\n');
fprintf('--    -------------    ----------    --------\n');

for r = ranks
    A_r = U(:, 1:r) * S(1:r, 1:r) * V(:, 1:r)';
    
    error_fro = norm(A - A_r, 'fro');
    error_2 = norm(A - A_r, 2);
    relative_error = error_fro / norm(A, 'fro');
    
    fprintf('%2d    %13.6f    %10.6f    %8.4f\n', ...
            r, error_fro, error_2, relative_error);
end

% 验证Eckart-Young定理
fprintf('\nEckart-Young定理验证:\n');
for r = [2, 4, 6]
    theoretical_error_2 = singular_values(r+1);
    A_r = U(:, 1:r) * S(1:r, 1:r) * V(:, 1:r)';
    actual_error_2 = norm(A - A_r, 2);
    
    fprintf('秩%d: 理论误差 = %.6f, 实际误差 = %.6f, 差异 = %.2e\n', ...
            r, theoretical_error_2, actual_error_2, ...
            abs(theoretical_error_2 - actual_error_2));
end

%% 2. 随机化SVD
fprintf('\n2. 随机化SVD\n');
fprintf('目标：快速计算低秩近似\n\n');

% 大矩阵测试
m_large = 500; n_large = 300; target_rank = 20;
A_large = randn(m_large, target_rank) * randn(target_rank, n_large) + ...
          0.1 * randn(m_large, n_large);

fprintf('大矩阵: %d×%d，目标秩 = %d\n', m_large, n_large, target_rank);

% 标准SVD（只计算前target_rank个）
tic;
[U_std, S_std, V_std] = svds(A_large, target_rank);
time_std = toc;

% 随机化SVD实现
tic;
[U_rand, S_rand, V_rand] = randomized_svd(A_large, target_rank, target_rank + 10);
time_rand = toc;

fprintf('计算时间比较:\n');
fprintf('标准SVD: %.4f 秒\n', time_std);
fprintf('随机化SVD: %.4f 秒\n', time_rand);
fprintf('加速比: %.2f\n', time_std / time_rand);

% 精度比较
A_approx_std = U_std * S_std * V_std';
A_approx_rand = U_rand * S_rand * V_rand';

error_std = norm(A_large - A_approx_std, 'fro');
error_rand = norm(A_large - A_approx_rand, 'fro');

fprintf('近似精度比较:\n');
fprintf('标准SVD误差: %.6f\n', error_std);
fprintf('随机化SVD误差: %.6f\n', error_rand);
fprintf('相对精度损失: %.2f%%\n', 100 * (error_rand - error_std) / error_std);

%% 3. 增量SVD
fprintf('\n3. 增量SVD\n');
fprintf('目标：动态更新低秩近似\n\n');

% 初始矩阵
A_base = randn(100, 80);
k = 10;  % 保持的秩

fprintf('初始矩阵: %d×%d，保持秩 = %d\n', size(A_base), k);

% 初始SVD
[U_base, S_base, V_base] = svds(A_base, k);
fprintf('初始SVD完成\n');

% 添加新列
new_cols = randn(100, 5);
A_extended = [A_base, new_cols];

fprintf('添加 %d 列后: %d×%d\n', size(new_cols, 2), size(A_extended));

% 直接SVD
tic;
[U_direct, S_direct, V_direct] = svds(A_extended, k);
time_direct = toc;

% 增量更新
tic;
[U_incremental, S_incremental, V_incremental] = incremental_svd_update(...
    U_base, S_base, V_base, new_cols, k);
time_incremental = toc;

fprintf('计算时间比较:\n');
fprintf('直接SVD: %.4f 秒\n', time_direct);
fprintf('增量SVD: %.4f 秒\n', time_incremental);
fprintf('加速比: %.2f\n', time_direct / time_incremental);

% 精度比较
A_approx_direct = U_direct * S_direct * V_direct';
A_approx_incremental = U_incremental * S_incremental * V_incremental';

error_direct = norm(A_extended - A_approx_direct, 'fro');
error_incremental = norm(A_extended - A_approx_incremental, 'fro');

fprintf('近似精度比较:\n');
fprintf('直接SVD误差: %.6f\n', error_direct);
fprintf('增量SVD误差: %.6f\n', error_incremental);

%% 4. 非负矩阵分解 (NMF)
fprintf('\n4. 非负矩阵分解 (NMF)\n');
fprintf('目标：约束低秩近似\n\n');

% 构造非负矩阵
m_nmf = 50; n_nmf = 40; r_nmf = 8;
W_true = abs(randn(m_nmf, r_nmf));
H_true = abs(randn(r_nmf, n_nmf));
A_nmf = W_true * H_true + 0.1 * abs(randn(m_nmf, n_nmf));

fprintf('非负矩阵: %d×%d，目标秩 = %d\n', m_nmf, n_nmf, r_nmf);
fprintf('矩阵元素范围: [%.3f, %.3f]\n', min(A_nmf(:)), max(A_nmf(:)));

% NMF算法（乘性更新）
max_iter = 100;
tol = 1e-6;

% 初始化
W_nmf = abs(randn(m_nmf, r_nmf));
H_nmf = abs(randn(r_nmf, n_nmf));

fprintf('NMF迭代过程:\n');
fprintf('迭代    目标函数值    相对变化\n');
fprintf('----    ----------    --------\n');

obj_prev = norm(A_nmf - W_nmf * H_nmf, 'fro')^2;

for iter = 1:max_iter
    % 更新H
    H_nmf = H_nmf .* (W_nmf' * A_nmf) ./ (W_nmf' * W_nmf * H_nmf + eps);
    
    % 更新W
    W_nmf = W_nmf .* (A_nmf * H_nmf') ./ (W_nmf * H_nmf * H_nmf' + eps);
    
    % 计算目标函数
    obj_current = norm(A_nmf - W_nmf * H_nmf, 'fro')^2;
    relative_change = abs(obj_current - obj_prev) / obj_prev;
    
    if mod(iter, 10) == 0 || iter <= 5
        fprintf('%4d    %10.6f    %8.2e\n', iter, obj_current, relative_change);
    end
    
    if relative_change < tol
        fprintf('收敛于迭代 %d\n', iter);
        break;
    end
    
    obj_prev = obj_current;
end

% NMF结果分析
A_nmf_approx = W_nmf * H_nmf;
nmf_error = norm(A_nmf - A_nmf_approx, 'fro');

% 与SVD比较
[U_svd, S_svd, V_svd] = svds(A_nmf, r_nmf);
A_svd_approx = U_svd * S_svd * V_svd';
svd_error = norm(A_nmf - A_svd_approx, 'fro');

fprintf('近似质量比较:\n');
fprintf('NMF误差: %.6f\n', nmf_error);
fprintf('SVD误差: %.6f\n', svd_error);
fprintf('NMF相对误差增加: %.2f%%\n', 100 * (nmf_error - svd_error) / svd_error);

% 检查非负性
fprintf('非负性检查:\n');
fprintf('W最小值: %.6f\n', min(W_nmf(:)));
fprintf('H最小值: %.6f\n', min(H_nmf(:)));

%% 5. 稀疏低秩近似
fprintf('\n5. 稀疏低秩近似\n');
fprintf('目标：同时实现低秩和稀疏性\n\n');

% 构造稀疏+低秩矩阵
m_sparse = 60; n_sparse = 50;
rank_sparse = 5;
sparsity_level = 0.1;

% 低秩部分
L_true = randn(m_sparse, rank_sparse) * randn(rank_sparse, n_sparse);

% 稀疏部分
S_true = zeros(m_sparse, n_sparse);
sparse_indices = randperm(m_sparse * n_sparse, round(sparsity_level * m_sparse * n_sparse));
S_true(sparse_indices) = 5 * randn(length(sparse_indices), 1);

% 观测矩阵
A_sparse = L_true + S_true;

fprintf('稀疏+低秩矩阵: %d×%d\n', m_sparse, n_sparse);
fprintf('真实秩: %d\n', rank_sparse);
fprintf('稀疏度: %.1f%%\n', 100 * sparsity_level);

% 简化的稀疏低秩分解（交替最小化）
max_iter_sparse = 50;
lambda = 0.1;  % 稀疏正则化参数

% 初始化
L_est = zeros(m_sparse, n_sparse);
S_est = A_sparse;

fprintf('稀疏低秩分解迭代:\n');
fprintf('迭代    ||L||_*    ||S||_1    重构误差\n');
fprintf('----    -------    -------    --------\n');

for iter = 1:max_iter_sparse
    % 更新L（低秩部分）
    [U_temp, S_temp, V_temp] = svd(A_sparse - S_est);
    s_temp = diag(S_temp);
    s_temp = max(s_temp - lambda, 0);  % 软阈值
    r_temp = sum(s_temp > 0);
    if r_temp > 0
        L_est = U_temp(:, 1:r_temp) * diag(s_temp(1:r_temp)) * V_temp(:, 1:r_temp)';
    else
        L_est = zeros(m_sparse, n_sparse);
    end
    
    % 更新S（稀疏部分）
    residual = A_sparse - L_est;
    S_est = sign(residual) .* max(abs(residual) - lambda, 0);  % 软阈值
    
    % 计算指标
    nuclear_norm = sum(s_temp);
    l1_norm = sum(abs(S_est(:)));
    reconstruction_error = norm(A_sparse - L_est - S_est, 'fro');
    
    if mod(iter, 10) == 0 || iter <= 5
        fprintf('%4d    %7.3f    %7.3f    %8.4f\n', ...
                iter, nuclear_norm, l1_norm, reconstruction_error);
    end
    
    if reconstruction_error < 1e-6
        fprintf('收敛于迭代 %d\n', iter);
        break;
    end
end

% 分解质量评估
L_error = norm(L_true - L_est, 'fro') / norm(L_true, 'fro');
S_error = norm(S_true - S_est, 'fro') / norm(S_true, 'fro');

fprintf('分解质量:\n');
fprintf('低秩部分相对误差: %.4f\n', L_error);
fprintf('稀疏部分相对误差: %.4f\n', S_error);

%% 6. 矩阵补全
fprintf('\n6. 矩阵补全\n');
fprintf('目标：从部分观测恢复低秩矩阵\n\n');

% 构造低秩矩阵
m_mc = 40; n_mc = 30; rank_mc = 6;
U_mc_true = randn(m_mc, rank_mc);
V_mc_true = randn(n_mc, rank_mc);
M_true = U_mc_true * V_mc_true';

% 随机观测
observation_rate = 0.3;
Omega = rand(m_mc, n_mc) < observation_rate;
M_observed = M_true;
M_observed(~Omega) = 0;

fprintf('矩阵补全问题: %d×%d，真实秩 = %d\n', m_mc, n_mc, rank_mc);
fprintf('观测率: %.1f%% (%d/%d 元素)\n', ...
        100 * observation_rate, sum(Omega(:)), numel(Omega));

% 简化的矩阵补全算法（交替最小二乘）
max_iter_mc = 100;
rank_estimate = 10;  % 估计秩

% 初始化
U_mc = randn(m_mc, rank_estimate);
V_mc = randn(n_mc, rank_estimate);

fprintf('矩阵补全迭代:\n');
fprintf('迭代    观测误差    总体误差    相对变化\n');
fprintf('----    --------    --------    --------\n');

obj_prev_mc = inf;

for iter = 1:max_iter_mc
    % 固定V，更新U
    for i = 1:m_mc
        observed_cols = find(Omega(i, :));
        if ~isempty(observed_cols)
            V_obs = V_mc(observed_cols, :);
            m_obs = M_observed(i, observed_cols)';
            U_mc(i, :) = (V_obs' * V_obs) \ (V_obs' * m_obs);
        end
    end
    
    % 固定U，更新V
    for j = 1:n_mc
        observed_rows = find(Omega(:, j));
        if ~isempty(observed_rows)
            U_obs = U_mc(observed_rows, :);
            m_obs = M_observed(observed_rows, j);
            V_mc(j, :) = (U_obs' * U_obs) \ (U_obs' * m_obs);
        end
    end
    
    % 计算误差
    M_reconstructed = U_mc * V_mc';
    observed_error = norm((M_true - M_reconstructed) .* Omega, 'fro') / ...
                     norm(M_true .* Omega, 'fro');
    total_error = norm(M_true - M_reconstructed, 'fro') / norm(M_true, 'fro');
    
    obj_current = norm((M_observed - M_reconstructed) .* Omega, 'fro')^2;
    relative_change = abs(obj_current - obj_prev_mc) / abs(obj_prev_mc);
    
    if mod(iter, 20) == 0 || iter <= 5
        fprintf('%4d    %8.4f    %8.4f    %8.2e\n', ...
                iter, observed_error, total_error, relative_change);
    end
    
    if relative_change < 1e-6
        fprintf('收敛于迭代 %d\n', iter);
        break;
    end
    
    obj_prev_mc = obj_current;
end

% 补全质量评估
final_error = norm(M_true - M_reconstructed, 'fro') / norm(M_true, 'fro');
fprintf('最终补全相对误差: %.4f\n', final_error);

% 估计秩分析
[~, S_mc, ~] = svd(M_reconstructed);
estimated_singular_values = diag(S_mc);
fprintf('重构矩阵的前8个奇异值: [%.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f, %.3f]\n', ...
        estimated_singular_values(1:8));

%% 总结
fprintf('\n=== 低秩近似算法总结 ===\n');
fprintf('1. SVD最优近似：理论最优，计算复杂度高\n');
fprintf('2. 随机化SVD：快速近似，适合大规模问题\n');
fprintf('3. 增量SVD：动态更新，适合流数据\n');
fprintf('4. 非负矩阵分解：约束近似，解释性好\n');
fprintf('5. 稀疏低秩分解：同时建模稀疏和低秩结构\n');
fprintf('6. 矩阵补全：从部分观测恢复完整矩阵\n\n');

fprintf('算法选择指南：\n');
fprintf('• 精度优先：标准SVD\n');
fprintf('• 速度优先：随机化SVD\n');
fprintf('• 非负约束：NMF\n');
fprintf('• 稀疏结构：稀疏低秩分解\n');
fprintf('• 缺失数据：矩阵补全\n');
fprintf('• 动态数据：增量SVD\n\n');

end

% 辅助函数
function [U, S, V] = randomized_svd(A, k, p)
% 随机化SVD实现
if nargin < 3
    p = min(k + 5, size(A, 2));
end

% 随机投影
Omega = randn(size(A, 2), p);
Y = A * Omega;

% QR分解
[Q, ~] = qr(Y, 0);

% 投影
B = Q' * A;

% SVD分解
[U_tilde, S, V] = svd(B, 'econ');

% 恢复U
U = Q * U_tilde;

% 截断到k个成分
if k < size(S, 1)
    U = U(:, 1:k);
    S = S(1:k, 1:k);
    V = V(:, 1:k);
end
end

function [U_new, S_new, V_new] = incremental_svd_update(U, S, V, new_cols, k)
% 增量SVD更新（简化版本）
m = size(U, 1);
r = size(S, 1);
p = size(new_cols, 2);

% 投影新列到当前子空间
proj_cols = U' * new_cols;
residual_cols = new_cols - U * proj_cols;

% QR分解残差
[Q_res, R_res] = qr(residual_cols, 0);

% 构造扩展矩阵
S_extended = [S, proj_cols; zeros(size(Q_res, 2), r), R_res];

% SVD分解
[U_tilde, S_new, V_tilde] = svd(S_extended);

% 更新U和V
U_extended = [U, Q_res];
U_new = U_extended * U_tilde;

V_extended = [V, zeros(size(V, 1), p); zeros(p, r), eye(p)];
V_new = V_extended * V_tilde;

% 截断到k个成分
if k < size(S_new, 1)
    U_new = U_new(:, 1:k);
    S_new = S_new(1:k, 1:k);
    V_new = V_new(:, 1:k);
end
end